knitr::opts_chunk$set(
echo = TRUE,
error = FALSE,
message = FALSE,
warning = FALSE
)
library(afex)
library(tidyverse)
library(kableExtra)
library(knitr)
library(plotly)
library(broom)
library(wesanderson)
library(robust)
library(ggparl)
library(mediation)
set.seed(2022)# reads excel data sheet
data <- read.csv("../data/study4.csv") # reads data
data_sat <- read.csv("../data/raw/stage3_satisfaction.csv") # reads in satisfaction data
# summary(data) # explore the data types, missing data and typos
data <- data %>%
mutate(gender = ifelse(gender == "m", "m", "f")) # corrects an issuewith extra space in "f "
data$gender <- as.factor(data$gender)
data$pp_number <- as.factor(data$pp_number)
data$condition <-as.factor(data$condition)
# mood data
mood1 <- read.csv("../data/raw/stage1_mood.csv") %>%
dplyr::select(subject, stimulusitem2, response)%>%
rename(response1=response, feeling = stimulusitem2)%>%
mutate(stage = 1)%>%
filter(subject != 6666)
mood2 <- read.csv("../data/raw/stage2_mood.csv") %>%
dplyr::select(subject, stimulusitem2, response)%>%
rename(response2=response, feeling = stimulusitem2)%>%
mutate(stage = 2)%>%
filter(subject != 6666)
mood3 <- read.csv("../data/raw/stage3_mood.csv") %>%
dplyr::select(subject, stimulusitem2, response) %>%
rename(response3=response, feeling = stimulusitem2)%>%
mutate(stage = 3)%>%
filter(subject != 6666)
mood_dat <- inner_join( mood1, mood2,
by = c("subject", "feeling"),
keep = FALSE)%>%
inner_join(., mood3,
by = c("subject", "feeling"),
keep = FALSE)%>%
dplyr::select(- starts_with("stage")) %>%
pivot_longer(cols = response1:response3,
names_to ="stage",
values_to = "rating")%>%
pivot_wider(names_from = "feeling",
values_from ="rating" )%>%
mutate(alert = (alert +strong + `clear-headed` + `well-coordinated` + energetic + `quick-witted` + attentive + proficient + interested )/9,
content =(contented + tranquil + amicable + gregarious)/4,
calm =(calm + relaxed)/2)
#`clear-headed` + `well-coordinated` + energetic + `quick-witted` + attentive + proficient + interested
# qualtrics data
#mainly info about drinking
data_q<- read.csv("../data/other/study_4Q.csv")%>%
dplyr::select(Q21, gender, age, Q15, Q16)%>%
rename(subject = Q21,
units = Q15,
units_beer = Q16)%>%
mutate(subject = as.numeric(subject),
condition = subject %/% 1000)%>%
filter(!is.na(subject))%>%
mutate(subject = as.factor(subject),
gender=as.factor(gender),
age = as.numeric(age),
units = as.numeric(units),
units_beer = as.numeric(units_beer))%>%
filter( subject != "6666")
# expectations & perception of taste/flavour
taste_dat <- read.csv("../data/raw/stage1_perception.csv")%>%
dplyr::select(subject, trialcode, response)%>%
rename(Epercept = trialcode, expected = response)
taste_dat$Epercept<- as.factor(taste_dat$Epercept)
expect_dat<- read.csv("../data/raw/stage1_expectations.csv")%>%
dplyr::select(subject, trialcode, response, stimulusitem3)%>%
filter(stimulusitem3 == "extremely")%>%
dplyr::select(- stimulusitem3)%>%
rename(percept = trialcode, perceived = response)
expect_dat$percept<- as.factor(expect_dat$percept)
# satisfaction data
#merge data_sat and expect_dat
satis_datP <- data_sat %>%
dplyr::select(subject, trialcode, response)%>%
dplyr::filter(trialcode == "VAS")
satis_datE <- expect_dat %>%
dplyr::filter(percept == "VAS_satisE_VAS")
wtp_dat <- data_sat %>%
filter(trialcode == "VAS_wtp")%>%
dplyr::select(subject, response)
# ITT data
itt_s1 <-read.csv("../data/summaries/stage1_ITT.csv")%>%
dplyr::select(script.subjectid, expressions.propcorrect_stim1:expressions.propcorrect_stim13)%>%
mutate(stage1 = 1)%>%
rename(subject = script.subjectid)%>%
rename_at(vars(matches("^expressions")), ~ str_remove(.,"^expressions.propcorrect_")) %>%
pivot_longer(cols = starts_with("stim") ,
names_to = "stimulus" ,
values_to = "correct")
itt_s2 <-read.csv("../data/summaries/stage2_ITT.csv")%>%
dplyr::select(script.subjectid, expressions.propcorrect_stim1:expressions.propcorrect_stim13)%>%
mutate(stage2 = 2)%>%
rename(subject = script.subjectid)%>%
rename_at(vars(matches("^expressions")), ~ str_remove(.,"^expressions.propcorrect_"))%>%
pivot_longer(cols = starts_with("stim") ,
names_to = "stimulus" ,
values_to = "correct")
itt_s3 <-read.csv("../data/summaries/stage3_ITT.csv")%>%
dplyr::select(script.subjectid, expressions.propcorrect_stim1:expressions.propcorrect_stim13)%>%
mutate(stage3 = 3)%>%
rename(subject = script.subjectid)%>%
rename_at(vars(matches("^expressions")), ~ str_remove(.,"^expressions.propcorrect_"))%>%
pivot_longer(cols = starts_with("stim") ,
names_to = "stimulus" ,
values_to = "correct")
itt_dat <- inner_join(itt_s1, itt_s2,
by = c("subject","stimulus"),
keep = FALSE,
suffix = c(".1", ".2"))%>%
inner_join(., itt_s3,
by = c("subject", "stimulus"),
keep = FALSE) %>%
pivot_longer(cols = starts_with("correct") ,
names_to = "stage",
values_to = "correct")%>%
dplyr::select(-c( stage1, stage2, stage3))%>%
mutate( stage = case_when(stage == "correct.1" ~ 'stage 1',
stage == "correct.2" ~ "stage 2",
stage == "correct" ~ "stage 3"))%>%
mutate(condition = subject %/% 1000) %>%
mutate(alcohol = if_else((condition%%2)==1, "alcoholic", "non-alcoholic"))%>%
mutate(focus = case_when(condition == 1 ~ 'control',
condition == 2 ~ 'control',
condition == 3 ~ 'hedonic',
condition == 4 ~ 'hedonic',
condition == 5 ~ 'utilitarian',
condition == 6 ~ 'utilitarian'))DA <- data %>%
dplyr::select(pp_number, BMI)%>%
rename(subject = pp_number)
DB <- data_sat %>%
dplyr::select(subject, response, trialcode ) %>%
mutate(row = row_number(), subject = as.factor(subject)) %>%
pivot_wider(names_from = trialcode, values_from = response) %>%
rename( wtp = VAS_wtp)%>%
dplyr::select(- row)
D1 <- inner_join( DA,DB,
by = "subject",
keep = FALSE
) %>%
pivot_longer(cols =c(VAS, wtp) ,
names_to = "measure",
values_to = "rating") %>%
inner_join(., data_q,
by = "subject",
keep = FALSE)%>%
filter(!is.na(rating))
# long dataset incl,. info about satisfaction and willingness to pay by condition, and alcohol/ non-alcohol
# need to add focus condition
#Joining mood data
D2 <- inner_join(mood1, mood2,
by = c("subject", "feeling"),
keep = FALSE) %>%
inner_join(., mood3,
by = c("subject","feeling"),
keep = FALSE)%>%
dplyr::select(- starts_with("stage")) %>%
mutate(condition = subject %/% 1000) %>%
pivot_longer(cols = starts_with("response"),
names_to = "stage",
names_prefix = "response",
values_to = "rating")
D3_dat <- inner_join(expect_dat, taste_dat,
by = "subject",
keep = FALSE)%>%
mutate(condition = subject %/% 1000) %>%
pivot_wider(names_from = percept,
values_from = perceived) %>%
pivot_wider(names_from = Epercept,
values_from = expected)%>%
rename_at(vars(matches("^VAS\\_")), ~ str_remove(.,"^VAS\\_"))%>%
rename_at(vars(matches("_VAS")), ~ str_remove(.,"_VAS")) %>%
mutate(alcohol = if_else((condition%%2)==1, "alcoholic", "non-alcoholic"))%>%
mutate(focus = case_when(condition == 1 ~ 'control',
condition == 2 ~ 'control',
condition == 3 ~ 'hedonic',
condition == 4 ~ 'hedonic',
condition == 5 ~ 'utilitarian',
condition == 6 ~ 'utilitarian'))
D3_dat$alcohol<- as.factor(D3_dat$alcohol)
D3_dat$focus<- as.factor(D3_dat$focus)
D3_dat$condition <- as.factor(D3_dat$condition)
satis_dat <- inner_join(satis_datE, satis_datP,
by = "subject",
keep = FALSE)%>%
dplyr::select(subject, perceived, response,)%>%
rename(satisE = perceived, satis = response)%>%
inner_join(., wtp_dat,
by = "subject",
keep = FALSE)%>%
rename(wtp = response) %>%
mutate(condition = subject %/% 1000) %>%
mutate(alcohol = if_else((condition%%2)==1, "alcoholic", "non-alcoholic"))%>%
mutate(focus = case_when(condition == 1 ~ 'control',
condition == 2 ~ 'control',
condition == 3 ~ 'hedonic',
condition == 4 ~ 'hedonic',
condition == 5 ~ 'utilitarian',
condition == 6 ~ 'utilitarian'))# gender
total_n <- nrow(data)
males <- data_q%>%filter(gender == "male") %>% nrow()
females <- data_q%>%filter(gender == "female") %>% nrow()
nonbinary <- data_q %>% filter(gender== "non-binary" )%>% nrow()
# by condition
gender_tib <- D1 %>%
filter(measure == "wtp") %>% # avoids duplicating data
group_by( condition) %>%
summarise(n=n(),
BMI = mean(BMI, na.rm = TRUE),
age = mean(age, na.rm = TRUE)) %>%
kable()xtabs(~ gender + condition , data = data_q,drop.unused.levels = TRUE)%>%
kable(digits = 2, caption = "A contingency table with observed gender values/counts across conditions") %>%
kable_styling(full_width = TRUE,
position = "left")| 1 | 2 | 3 | 4 | 5 | 6 | |
|---|---|---|---|---|---|---|
| female | 24 | 21 | 22 | 25 | 30 | 23 |
| male | 6 | 8 | 7 | 6 | 6 | 8 |
| non-binary | 2 | 1 | 1 | 1 | 0 | 0 |
chi_test <- chisq.test(data_q$gender, data_q$condition, simulate.p.value = TRUE)
chi_test$expected %>%
kable(digits = 1, caption = "A contingency table of expected gender values/counts across conditions")%>%
kable_styling(full_width = TRUE,
position = "left")| 1 | 2 | 3 | 4 | 5 | 6 | |
|---|---|---|---|---|---|---|
| female | 24.3 | 22.8 | 22.8 | 24.3 | 27.3 | 23.5 |
| male | 6.9 | 6.4 | 6.4 | 6.9 | 7.7 | 6.7 |
| non-binary | 0.8 | 0.8 | 0.8 | 0.8 | 0.9 | 0.8 |
# χ2 test requires all expected frequencies to be greater than five --> simulate.p.value avoids the problemAltogether 189 participants took part in the study, as is common for psychology studies, most of these were women (145 self-identified as female, 41 as male and 5 as non-binary). The gender split across the conditions was not problematic \(\chi^2\) = 5.28, p = 0.903 (simulated p values based on 2000 replicates, as expected values << 5).
# plot
## bar plot is best here
gender_pal <- c("#d62828", "#003049", "#9d69a3")
gender_cond_plot <-
ggplot(data = data_q)+
geom_bar(aes(x= condition, fill = gender))+
scale_fill_manual(values = gender_pal) +
scale_y_continuous(name="", limits=c(0, 40), breaks = seq(0,40, 5)) +
scale_x_discrete(name = "conditions", breaks = seq(1,6,1), limits = c("1", "2", "3", "4", "5", "6"))+
theme(panel.grid = element_blank())+
theme_minimal()
ggplotly(gender_cond_plot) Number of participants in each conditions
data_desc <- D1 %>%
filter(measure == "wtp")
age_lm <- lm(age ~ condition, data=data_desc)
age_sum <- broom::tidy(age_lm)
age_glance <- broom::glance(age_lm)Participants’ age ranged between 18 and 57. Most participants were very young and the mean age of participants was 20.3 (sd = 4.7). The age did not significantly differ across the six conditions F(181, 1) = 0.26, p = 0.608.
age_plot <- ggplot(data= data_desc, aes(x=condition, y = age, group = condition) )+
geom_boxplot(fill= "#6c757d", alpha = 0.5, outlier.shape = NA) +
geom_jitter(shape=16, position=position_jitter(0.2), aes(color = gender))+
scale_y_continuous(name="age", limits=c(18, 60), breaks = seq(0,60, 10)) +
scale_x_discrete(name = "conditions", limits = c("1", "2", "3", "4", "5", "6"))+
scale_color_manual(values = gender_pal) +
theme_minimal()
age_plot2 <- ggplot(data= data_desc, aes(x=condition, y = age, group = condition) )+
geom_boxplot(fill= "#6c757d", alpha = 0.5, outlier.shape = NA) +
geom_jitter(shape=16, position=position_jitter(0.2), aes(color = gender))+
scale_y_continuous(name="age", limits=c(15, 30), breaks = seq(15,30,5)) +
scale_x_discrete(name = "conditions", limits = c("1", "2", "3", "4", "5", "6"))+
scale_color_manual(values = gender_pal) +
theme_minimal()
ggplotly(age_plot)Participants’ age by condition and gender. Note range of the y-axis
ggplotly(age_plot2)Participants’ age by condition and gender. Note range of the y-axis
bmi_lm <- lm(BMI ~ condition, data=data_desc)
bmi_glance <- broom::glance(bmi_lm)The participants’ BMI ranged between 16.7 - 36.1. The average BMI was 23.1 (sd =3.7) and again, this did not differ across conditions F(186, 1) = 1.62, p = 0.205.
BMI_plot <- ggplot(aes(x=condition, y = BMI), data= data_desc)+
geom_boxplot(fill= "#6c757d", alpha = 0.5) +
geom_jitter(shape=16, position=position_jitter(0.2), aes(color = gender))+
scale_y_continuous(name="BMI", limits=c(18, 35), breaks = seq(0,35, 10)) +
scale_x_discrete(name = "conditions", limits = c("1", "2", "3", "4", "5", "6"))+
scale_color_manual(values = gender_pal) +
theme_minimal()
ggplotly(BMI_plot)Participants’ BMI by condition and gender
Visual inspection of participants’ baseline mood does not give raise to any concerns. Do I need to run 16 anovas? Does not seem necessary…)
mood_base <- D2 %>%
filter(stage == 1)
palette <- c("#033f63","#28666e","#7c9885","#b5b682","#fedc97","#e5e1ee","#c17767","#210203","#eb9486","#82204a", "#0a2342","#2ca58d","#84bc9c","#fffdf7","#f46197","#fec601")
ggplot(data = mood_base, aes(x=condition, y = rating, group = condition, fill = feeling))+
geom_boxplot()+
facet_wrap("feeling")+
scale_fill_manual(values = palette) +
scale_y_continuous(name="", limits=c(0, 100), breaks = seq(0,100, 20)) +
theme(panel.grid = element_blank())+
theme_minimal()
Participants’ baseline mood ratings
While the average alcohol consumption was realtively low 10.6 (sd = 8.9), there was a considerable differences among participants. The lowest alcohol consumption was only 1, equivalent to a half a pint of light beer (3.5%) once a week, the highest alcohol consumption was 70 units of alcohol a week, which is equivalent to over 20 pints of strong (~5% abv) beer every week. Interestingly, 45 (23.8%) participants consumed more than the recommended 14 units a week. We also asked participants to estimate how many units of beer they consumed. This was on average 4.8 (sd = 4.7), this however varied between 0 and 30. Thefigures below show participants’ weekly alcohol consumption, weekly beer consumption and the proportion of units of alcohol that came from beer. Neither the overall alcohol consumption nor beer consumption differed significantly across the conditions, F(2, 189) = 3.19, p = 0.069 and F(2, 189) = 0.06, p = 0.804, respectively.
# plot units alcohol consumption by condition
units_plot<-ggplot(aes(x=condition,
y = units,
group = condition),
data= data_q)+
geom_hline(aes(yintercept = 14, linetype = "14 units/week"),
color = "red")+
scale_linetype_manual(name = "max. recommended\nalcohol consumption", values = c(2, 2),
guide = guide_legend(override.aes = list(color = c("red")), title.theme = element_text(size = 8, face = "bold" ) ))+
geom_boxplot(fill= "#6c757d",
alpha = 0.5,
outlier.shape = NA) +
geom_jitter(shape=16,
position=position_jitter(0.2),
aes(color = gender))+
scale_y_continuous(name="alcohol consumption\n(units)",
limits=c(0, 75),
breaks = seq(0,75, 5)) +
scale_x_discrete(name = "conditions",
limits = c("1", "2", "3", "4", "5", "6"))+
scale_color_manual(values = gender_pal) +
theme_minimal()
ggplotly(units_plot)participants’ alcohol consumption by condition
# plot units beer consumption by condition
beer_plot<-ggplot(aes(x=condition,
y = units_beer,
group = condition),
data= data_q)+
geom_boxplot(fill= "#6c757d",
alpha = 0.5,
outlier.shape = NA) +
geom_jitter(shape=16,
position=position_jitter(0.2),
aes(color = gender))+
scale_y_continuous(name="beer consumption\n(units)",
limits=c(0, 35),
breaks = seq(0,35, 5)) +
scale_x_discrete(name = "conditions",
limits = c("1", "2", "3", "4", "5", "6"))+
scale_color_manual(values = gender_pal) +
theme_minimal()
ggplotly(beer_plot)participants’ beer consumption by condition
summary_beer <- data_q %>%
pivot_longer(cols = starts_with("units"), names_to = "al_id", values_to = "units") %>%
mutate(al_id = ifelse(al_id =="units", "all", "beer"))%>%
group_by(condition, al_id)%>%
summarize(units = mean(units))%>%
pivot_wider(names_from = al_id, values_from = units)%>%
mutate(perc= (beer/all)*100)%>%
pivot_longer(cols = c(all, beer), names_to = "alcohol", values_to = "units")
al_pal <- c( "#1d3557","#f4a261")
beer_prop <- ggplot(data=summary_beer,
aes(x=condition,
y=units,
fill = alcohol)) +
geom_bar(stat="identity",
position=position_dodge())+
geom_text(aes(label =paste0(round(perc,1),"%") ),
data = summary_beer %>%filter(alcohol == "beer"),
show.legend = FALSE,
size = 2,
nudge_x = 0.3,
nudge_y = 1,
color = "#1d3557" )+
scale_fill_manual(values = al_pal)+
scale_y_continuous(name="alcohol consumption\n(units)",
limits=c(0, 15),
breaks = seq(0,15, 5)) +
scale_x_discrete(name = "conditions",
limits = c("1", "2", "3", "4", "5", "6"))+
labs( fill = "alcohol")+
theme_minimal()
ggplotly(beer_prop)% of alcohol consumption that is beer by condition
mediation package, function mediate() some resources and documentation
bitter_dat <- D3_dat%>%
dplyr::select(subject, alcohol, focus, condition, bitterE, bitter)
# prep dataset
# model m A --> M
mbm <- lm(bitterE ~ focus*alcohol, data = bitter_dat)
# model A--> B
mby <-lm(bitter ~ bitterE + focus*alcohol, data = bitter_dat)
# robustSE = TRUE to get heteroscedasticity consistent SE
# nonparametric bootstrap rather than the quasi-Bayesian Monte
# Carlo simulation for variance estimation via the boot = TRUE argument
# also entering alcohol/focus as a covariate
broom::tidy(mbm)## # A tibble: 6 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 55.5 2.80 19.8 1.61e-47
## 2 focushedonic 5.26 4.06 1.29 1.97e- 1
## 3 focusutilitarian -3.08 3.85 -0.801 4.24e- 1
## 4 alcoholnon-alcoholic -3.70 4.03 -0.919 3.59e- 1
## 5 focushedonic:alcoholnon-alcoholic 1.22 5.72 0.214 8.31e- 1
## 6 focusutilitarian:alcoholnon-alcoholic 6.99 5.59 1.25 2.13e- 1med_bitter_focusH <- mediate(mbm, mby,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "bitterE",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "hedonic",)
med_bitter_focusU <- mediate(mbm, mby,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "bitterE",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "utilitarian")
summary(med_bitter_focusH)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 1.5961 0.0275 3.82 0.044 *
## ADE 1.6530 -5.8790 9.24 0.678
## Total Effect 3.2490 -4.2910 10.74 0.422
## Prop. Mediated 0.4912 -3.3536 3.22 0.446
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_bitter_focusH)
summary(med_bitter_focusU)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.0925 -1.5001 1.73 0.88
## ADE 3.7717 -3.0102 10.89 0.32
## Total Effect 3.8642 -3.2772 10.99 0.29
## Prop. Mediated 0.0239 -1.5064 1.58 0.86
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_bitter_focusU)
med_bitter_alcohol <- mediate(mbm, mby,
sims = 1000,
boot.ci.type = "perc",
treat = "alcohol",
mediator = "bitterE",
robustSE = TRUE,
boot = TRUE,
conf.level = 0.95,
covariates = "focus",
control.value = "alcoholic",
treat.value = "non-alcoholic")
summary(med_bitter_alcohol)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.229 -1.629 1.14 0.70
## ADE -0.167 -6.500 6.16 0.99
## Total Effect -0.397 -6.912 5.91 0.92
## Prop. Mediated 0.578 -1.912 2.72 0.87
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_bitter_alcohol)
bitter_tib <-bitter_dat %>%
pivot_longer(cols = starts_with("bitter"),
names_to = "bitter",
values_to = "rating")%>%
group_by(alcohol, focus,bitter)%>%
summarise(n = n(),
meanB = mean(rating, na.rm = TRUE),
sdB = sd(rating, na.rm = TRUE),
ci_lo = meanB - 1.96 * sdB/sqrt(n),
ci_hi = meanB + 1.96 * sdB/sqrt(n))%>%
filter(!is.na(focus))
rating_labs <- c("perceived", "expected")
names(rating_labs) <- c("bitter", "bitterE")
bitter_plot<-ggplot(aes(x= focus, y = meanB, fill = alcohol), data = bitter_tib ) +
facet_grid(~ bitter,
labeller = labeller(bitter = rating_labs))+
geom_bar(stat="identity", position=position_dodge())+
geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=.2,
position=position_dodge(.9))+
scale_y_continuous(name="bitterness",
limits=c(0, 100),
breaks = seq(0,100, 25)) +
scale_x_discrete(name = "focus" )+
scale_fill_manual(values = gender_pal) +
theme_minimal()
ggplotly(bitter_plot)Expected and perceived ratings of bitterness errorbars indicate 95% CI (1.96 SE)
refresh_dat <- D3_dat%>%
dplyr::select(subject, alcohol, focus, condition, refreshmentE, refreshment)
# prep dataset
# model m A --> M
mrm <- lm(refreshmentE ~ focus*alcohol, data = refresh_dat)
# model A--> B
mry <-lm(refreshment ~ refreshmentE + focus*alcohol, data = refresh_dat)
broom::tidy(mrm)%>%
kable(digits = 3, caption = "effect of alcohol and focus on expected refreshment")%>%
kable_styling(full_width = TRUE)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 63.250 | 3.052 | 20.725 | 0.000 |
| focushedonic | 0.474 | 4.426 | 0.107 | 0.915 |
| focusutilitarian | 6.472 | 4.194 | 1.543 | 0.125 |
| alcoholnon-alcoholic | -2.417 | 4.387 | -0.551 | 0.582 |
| focushedonic:alcoholnon-alcoholic | 2.818 | 6.232 | 0.452 | 0.652 |
| focusutilitarian:alcoholnon-alcoholic | -0.435 | 6.095 | -0.071 | 0.943 |
med_refreshment_focusH <- mediate(mrm, mry,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "refreshmentE",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "hedonic")
summary(med_refreshment_focusH)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.476 -1.180 1.99 0.518
## ADE -5.184 -10.226 0.31 0.064 .
## Total Effect -4.708 -9.785 0.80 0.090 .
## Prop. Mediated -0.101 -1.410 0.73 0.584
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_refreshment_focusH)
med_refreshment_focusU <- mediate(mrm, mry,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "refreshmentE",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "utilitarian")
summary(med_refreshment_focusU)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 1.6072 0.0126 3.59 0.046 *
## ADE -3.8968 -8.8457 1.17 0.120
## Total Effect -2.2895 -7.5390 3.10 0.418
## Prop. Mediated -0.7020 -11.5112 7.92 0.456
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_refreshment_focusU)
med_refreshment_alcohol <- mediate(mrm, mry,
sims = 1000,
boot.ci.type = "perc",
treat = "alcohol",
mediator = "refreshmentE",
robustSE = TRUE,
boot = TRUE,
conf.level = 0.95,
covariates = "focus",
control.value = "alcoholic",
treat.value = "non-alcoholic")
summary(med_refreshment_alcohol)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.428 -1.945 0.80 0.50
## ADE -1.600 -5.955 2.61 0.43
## Total Effect -2.028 -6.515 2.32 0.33
## Prop. Mediated 0.211 -2.628 2.81 0.56
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_refreshment_alcohol)
refreshment_tib <-refresh_dat %>%
pivot_longer(cols = starts_with("refreshment"),
names_to = "refreshment",
values_to = "rating")%>%
group_by(alcohol, focus,refreshment)%>%
summarise(n = n(),
meanB = mean(rating, na.rm = TRUE),
sdB = sd(rating, na.rm = TRUE),
ci_lo = meanB - 1.96 * sdB/sqrt(n),
ci_hi = meanB + 1.96 * sdB/sqrt(n))%>%
filter(!is.na(focus))
rating_labs <- c("perceived", "expected")
names(rating_labs) <- c("refreshment", "refreshmentE")
refreshment_plot <- ggplot(aes(x= focus, y = meanB, fill = alcohol), data = refreshment_tib ) +
facet_grid(~ refreshment,
labeller = labeller(refreshment = rating_labs))+
geom_bar(stat="identity", position=position_dodge())+
geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=.2,
position=position_dodge(.9))+
scale_y_continuous(name="refreshment",
limits=c(0, 100),
breaks = seq(0,100, 25)) +
scale_x_discrete(name = "focus" )+
scale_fill_manual(values = gender_pal) +
theme_minimal()
ggplotly(refreshment_plot)Expected and perceived ratings of Refreshment errorbars indicate 95% CI (1.96 SE)
liking_dat <- D3_dat%>%
dplyr::select(subject, alcohol, focus, condition, likingE, liking)
# prep dataset
# model m A --> M
mlm <- lm(likingE ~ focus*alcohol, data = liking_dat)
# model A--> B
mly <-lm(liking ~ likingE + focus*alcohol, data = liking_dat)
broom::tidy(mlm)%>%
kable(digits = 3, caption = "effect of alcohol and focus on expected liking")%>%
kable_styling(full_width = TRUE)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 58.594 | 3.404 | 17.212 | 0.000 |
| focushedonic | -6.249 | 4.937 | -1.266 | 0.207 |
| focusutilitarian | -1.455 | 4.679 | -0.311 | 0.756 |
| alcoholnon-alcoholic | -5.394 | 4.894 | -1.102 | 0.272 |
| focushedonic:alcoholnon-alcoholic | 2.955 | 6.952 | 0.425 | 0.671 |
| focusutilitarian:alcoholnon-alcoholic | -0.261 | 6.798 | -0.038 | 0.969 |
med_liking_focusH <- mediate(mlm, mly,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "likingE",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "hedonic")
summary(med_liking_focusH)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.745 -2.151 0.39 0.202
## ADE -5.091 -11.422 1.42 0.112
## Total Effect -5.835 -12.054 0.77 0.078 .
## Prop. Mediated 0.128 -0.465 0.89 0.256
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_liking_focusH)
med_liking_focusU <- mediate(mlm, mly,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "likingE",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "utilitarian")
summary(med_liking_focusU)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.2454 -1.3588 0.92 0.67
## ADE -3.0973 -9.3443 2.94 0.30
## Total Effect -3.3427 -9.8839 2.99 0.27
## Prop. Mediated 0.0734 -1.3463 1.02 0.65
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_liking_focusU)
med_liking_alcohol <- mediate(mlm, mly,
sims = 1000,
boot.ci.type = "perc",
treat = "alcohol",
mediator = "likingE",
robustSE = TRUE,
boot = TRUE,
conf.level = 0.95,
covariates = "focus",
control.value = "alcoholic",
treat.value = "non-alcoholic")
summary(med_liking_alcohol)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.703 -2.170 0.11 0.12
## ADE -3.393 -8.464 2.01 0.20
## Total Effect -4.096 -9.395 1.24 0.11
## Prop. Mediated 0.172 -0.529 1.51 0.20
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_liking_alcohol)
liking_tib <-liking_dat %>%
pivot_longer(cols = starts_with("liking"),
names_to = "liking",
values_to = "rating")%>%
group_by(alcohol, focus,liking)%>%
summarise(n = n(),
meanB = mean(rating, na.rm = TRUE),
sdB = sd(rating, na.rm = TRUE),
ci_lo = meanB - 1.96 * sdB/sqrt(n),
ci_hi = meanB + 1.96 * sdB/sqrt(n))%>%
filter(!is.na(focus))
rating_labs <- c("perceived", "expected")
names(rating_labs) <- c("liking", "likingE")
liking_plot <- ggplot(aes(x= focus, y = meanB, fill = alcohol), data = liking_tib ) +
facet_grid(~ liking,
labeller = labeller(liking = rating_labs))+
geom_bar(stat="identity", position=position_dodge())+
geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=.2,
position=position_dodge(.9))+
scale_y_continuous(name="liking",
limits=c(0, 100),
breaks = seq(0,100, 25)) +
scale_x_discrete(name = "focus" )+
scale_fill_manual(values = gender_pal) +
theme_minimal()
ggplotly(liking_plot)Expected and perceived ratings of liking errorbars indicate 95% CI (1.96 SE)
body_dat <- D3_dat%>%
dplyr::select(subject, alcohol, focus, condition, bodyE, body)
# prep dataset
# model m A --> M
mbom <- lm(bodyE ~ focus*alcohol, data = body_dat)
# model A--> B
mboy <-lm(body ~ bodyE + focus*alcohol, data = body_dat)
broom::tidy(mbom)%>%
kable(digits = 3, caption = "effect of alcohol and focus on expected body") %>%
kable_styling(full_width = TRUE)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 48.594 | 3.100 | 15.676 | 0.000 |
| focushedonic | 1.682 | 4.496 | 0.374 | 0.709 |
| focusutilitarian | 3.240 | 4.260 | 0.760 | 0.448 |
| alcoholnon-alcoholic | -7.727 | 4.456 | -1.734 | 0.085 |
| focushedonic:alcoholnon-alcoholic | 6.857 | 6.330 | 1.083 | 0.280 |
| focusutilitarian:alcoholnon-alcoholic | 3.636 | 6.190 | 0.587 | 0.558 |
med_body_focusH <- mediate(mbom, mboy,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "bodyE",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "hedonic")
summary(med_body_focusH)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.996 -0.152 3.36 0.14
## ADE 1.278 -5.479 7.96 0.67
## Total Effect 2.273 -4.116 9.01 0.48
## Prop. Mediated 0.438 -2.955 2.96 0.53
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_body_focusH)
med_body_focusU <- mediate(mbom, mboy,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "bodyE",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "utilitarian")
summary(med_body_focusU)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.9918 -0.0983 3.05 0.10
## ADE -1.6508 -8.3809 5.66 0.64
## Total Effect -0.6590 -7.5640 6.56 0.86
## Prop. Mediated -1.5049 -4.2345 4.24 0.91
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_body_focusU)
med_body_alcohol <- mediate(mbom, mboy,
sims = 1000,
boot.ci.type = "perc",
treat = "alcohol",
mediator = "bodyE",
robustSE = TRUE,
boot = TRUE,
conf.level = 0.95,
covariates = "focus",
control.value = "alcoholic",
treat.value = "non-alcoholic")
summary(med_body_alcohol)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.8385 -2.5845 0.19 0.148
## ADE -6.1981 -12.1954 -0.50 0.042 *
## Total Effect -7.0366 -13.0981 -1.26 0.016 *
## Prop. Mediated 0.1192 -0.0543 0.58 0.156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 190
##
##
## Simulations: 1000plot(med_body_alcohol)
body_tib <-body_dat %>%
pivot_longer(cols = starts_with("body"),
names_to = "body",
values_to = "rating")%>%
group_by(alcohol, focus,body)%>%
summarise(n = n(),
meanB = mean(rating, na.rm = TRUE),
sdB = sd(rating, na.rm = TRUE),
ci_lo = meanB - 1.96 * sdB/sqrt(n),
ci_hi = meanB + 1.96 * sdB/sqrt(n))%>%
filter(!is.na(focus))
rating_labs <- c("perceived", "expected")
names(rating_labs) <- c("body", "bodyE")
body_plot <-ggplot(aes(x= focus, y = meanB, fill = alcohol), data = body_tib ) +
facet_grid(~ body,
labeller = labeller(body = rating_labs))+
geom_bar(stat="identity", position=position_dodge())+
geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=.2,
position=position_dodge(.9))+
scale_y_continuous(name="body",
limits=c(0, 100),
breaks = seq(0,100, 25)) +
scale_x_discrete(name = "focus" )+
scale_fill_manual(values = gender_pal) +
theme_minimal()
ggplotly(body_plot)Expected and perceived ratings of body errorbars indicate 95% CI (1.96 SE)
satis_tib <-satis_dat %>%
pivot_longer(cols = starts_with("satis"),
names_to = "satis",
values_to = "rating")%>%
group_by(alcohol, focus,satis)%>%
summarise(n = n(),
meanB = mean(rating, na.rm = TRUE),
sdB = sd(rating, na.rm = TRUE),
ci_lo = meanB - 1.96 * sdB/sqrt(n),
ci_hi = meanB + 1.96 * sdB/sqrt(n))%>%
filter(!is.na(focus))
rating_labs <- c("perceived", "expected")
names(rating_labs) <- c("satis", "satisE")
satis_plot<-ggplot(aes(x= focus, y = meanB, fill = alcohol), data = satis_tib ) +
facet_grid(~ satis,
labeller = labeller(satis = rating_labs))+
geom_bar(stat="identity", position=position_dodge())+
geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=.2,
position=position_dodge(.9))+
scale_y_continuous(name="satisfaction",
limits=c(0, 100),
breaks = seq(0,100, 25)) +
scale_x_discrete(name = "focus" )+
scale_fill_manual(values = gender_pal) +
theme_minimal()
ggplotly(satis_plot)Expected and perceived ratings of satisfaction errorbars indicate 95% CI (+/- 1.96 SE)
# model m A --> M
msm <- lm(satisE ~ focus*alcohol, data = satis_dat)
# model A--> B
msy <-lm(satis ~ satisE + focus*alcohol, data = satis_dat)
broom::tidy(msm)%>%
kable(digits = 3, caption = "effect of alcoholand focus on expected satisfaction")%>%
kable_styling(full_width = TRUE)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 55.500 | 3.088 | 17.972 | 0.000 |
| focushedonic | -1.893 | 4.444 | -0.426 | 0.671 |
| focusutilitarian | -0.203 | 4.155 | -0.049 | 0.961 |
| alcoholnon-alcoholic | -8.500 | 4.367 | -1.946 | 0.053 |
| focushedonic:alcoholnon-alcoholic | 5.955 | 6.183 | 0.963 | 0.337 |
| focusutilitarian:alcoholnon-alcoholic | 7.803 | 6.028 | 1.294 | 0.197 |
med_satis_focusU <- mediate(msm, msy,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "satisE",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "utilitarian")
summary(med_satis_focusU)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.682 -0.338 2.20 0.22
## ADE -0.476 -6.678 5.83 0.88
## Total Effect 0.206 -6.527 6.45 0.94
## Prop. Mediated 3.312 -3.778 4.27 0.88
##
## Sample Size Used: 187
##
##
## Simulations: 1000plot(med_satis_focusU)
med_satis_focusH <- mediate(msm, msy,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "satisE",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "hedonic")
summary(med_satis_focusH)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.194 -1.087 1.61 0.78
## ADE -4.819 -12.295 2.51 0.22
## Total Effect -4.625 -12.059 2.96 0.26
## Prop. Mediated -0.042 -1.181 1.09 0.90
##
## Sample Size Used: 187
##
##
## Simulations: 1000plot(med_satis_focusH)
med_satis_alcohol <- mediate(msm, msy,
sims = 1000,
boot.ci.type = "perc",
treat = "alcohol",
mediator = "satisE",
robustSE = TRUE,
boot = TRUE,
conf.level = 0.95,
covariates = "focus",
control.value = "alcoholic",
treat.value = "non-alcoholic")
summary(med_satis_alcohol)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.7113 -2.1322 0.17 0.150
## ADE -7.0777 -12.2471 -1.40 0.012 *
## Total Effect -7.7890 -13.0529 -2.07 0.008 **
## Prop. Mediated 0.0913 -0.0279 0.40 0.154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 187
##
##
## Simulations: 1000plot(med_satis_alcohol)
# model m A --> M
mwm <- lm(satis ~ focus*alcohol, data = satis_dat)
# model A--> B
mwy <-lm(wtp ~ satis + focus*alcohol, data = satis_dat)
med_wtp_focusH <- mediate(mwm, mwy,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "satis",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "hedonic")
summary(med_wtp_focusH)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0943 -0.2411 0.05 0.21
## ADE -0.0658 -0.3619 0.20 0.63
## Total Effect -0.1601 -0.4736 0.15 0.31
## Prop. Mediated 0.5889 -3.6505 4.28 0.36
##
## Sample Size Used: 187
##
##
## Simulations: 1000plot(med_wtp_focusH)
med_wtp_focusU <- mediate(mwm, mwy,
sims = 1000,
boot = TRUE,
boot.ci.type = "perc",
treat = "focus",
mediator = "satis",
covariates = "alcohol",
robustSE = TRUE,
conf.level = 0.95,
control.value = "control",
treat.value = "utilitarian")
summary(med_wtp_focusU)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.0042 -0.1337 0.14 0.95
## ADE -0.1707 -0.4372 0.10 0.20
## Total Effect -0.1665 -0.4779 0.15 0.28
## Prop. Mediated -0.0252 -3.1410 3.37 0.85
##
## Sample Size Used: 187
##
##
## Simulations: 1000plot(med_wtp_focusU)
med_wtp_alcohol <- mediate(mwm, mwy,
sims = 1000,
boot.ci.type = "perc",
treat = "alcohol",
mediator = "satis",
robustSE = TRUE,
boot = TRUE,
conf.level = 0.95,
covariates = "focus",
control.value = "alcoholic",
treat.value = "non-alcoholic")
summary(med_wtp_alcohol)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.1588 -0.2755 -0.04 0.004 **
## ADE -0.0949 -0.3235 0.13 0.402
## Total Effect -0.2536 -0.5000 -0.01 0.036 *
## Prop. Mediated 0.6260 0.1091 2.85 0.036 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 187
##
##
## Simulations: 1000plot(med_wtp_alcohol)
wtp_model <- robust::lmRob(wtp ~ focus * alcohol, data = satis_dat)
summary(wtp_model)##
## Call:
## robust::lmRob(formula = wtp ~ focus * alcohol, data = satis_dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1655 -0.7333 0.1019 0.8345 2.1912
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.1655 0.1790 17.683 <2e-16 ***
## focushedonic -0.2674 0.2554 -1.047 0.2966
## focusutilitarian -0.3567 0.2403 -1.484 0.1394
## alcoholnon-alcoholic -0.4321 0.2510 -1.721 0.0869 .
## focushedonic:alcoholnon-alcoholic 0.1903 0.3539 0.538 0.5915
## focusutilitarian:alcoholnon-alcoholic 0.3233 0.3474 0.931 0.3532
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9126 on 181 degrees of freedom
## Multiple R-Squared: 0.03504
##
## Test for Bias:
## statistic p-value
## M-estimate 14.972 0.02048
## LS-estimate -6.361 1.00000wtp_tib <- satis_dat %>%
group_by(focus, alcohol) %>%
summarise(n = n(),
mean_wtp = mean(wtp, na.rm = TRUE),
sd_wtp = sd(wtp, na.rm = TRUE),
ci_lo = mean_wtp - 1.96 * (sd_wtp/sqrt(n)),
ci_high = mean_wtp + 1.96 *(sd_wtp/sqrt(n)))
wtp_tib %>%
kable(digits = 2, caption = "Participants' mean willingness to pay in different alcohol and focus conditions.") %>%
kable_styling(full_width = TRUE)| focus | alcohol | n | mean_wtp | sd_wtp | ci_lo | ci_high |
|---|---|---|---|---|---|---|
| control | alcoholic | 30 | 3.13 | 0.90 | 2.81 | 3.46 |
| control | non-alcoholic | 30 | 2.73 | 0.83 | 2.44 | 3.03 |
| hedonic | alcoholic | 28 | 2.89 | 0.88 | 2.57 | 3.22 |
| hedonic | non-alcoholic | 32 | 2.66 | 0.97 | 2.32 | 2.99 |
| utilitarian | alcoholic | 37 | 2.84 | 0.73 | 2.60 | 3.07 |
| utilitarian | non-alcoholic | 30 | 2.70 | 0.88 | 2.39 | 3.01 |
wtp_bar <-ggplot(aes(x=focus,
y = mean_wtp,
fill = alcohol),
data= wtp_tib)+
geom_bar(stat="identity", position=position_dodge())+
geom_errorbar(aes(ymin=ci_lo, ymax=ci_high), width=.2,
position=position_dodge(.9))+
scale_y_continuous(name="willingness to pay (£)",
limits=c(0, 5),
breaks = seq(0,5, 1)) +
scale_x_discrete(name = "focus" )+
scale_fill_manual(values = gender_pal) +
theme_minimal()
wtp_bar
Participants’ willingness to pay for alcoholic and non-alcoholic beer. Error bars indicate confidence intervals (+/- 1.96 SE)
Participants’ cognitive performance was assessed using the inspection time task (10.1016/j.neuroimage.2004.03.047) which measures participants’information processing. Participants completed the task just before, 30 minutes and 60 minutes after consumption.
-pp performance improves as duration of stimulus increases - effect of alcohol, stimulus duration and stage - interactions: alcoholstimulus stimulusstage alcohol*stimulus*stage - no effect of focus
# these are two methods of analysis, should do the same thing, buttakes too long to run
# m_itt <- lmerTest::lmer(correct ~ focus * alcohol *stage *stimulus + (stage*stimulus| subject), data=itt_dat, REML=F)
#
# summary(m_itt)
# contrasts
itt_afx <- afex::aov_4(correct ~ focus*alcohol*stage*stimulus + (stimulus*stage|subject), data = itt_dat)
summary(itt_afx)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 5831.3 1 38.734 178 26797.4015 < 2.2e-16 ***
## focus 0.6 2 38.734 178 1.4063 0.247768
## alcohol 1.8 1 38.734 178 8.1225 0.004888 **
## focus:alcohol 0.7 2 38.734 178 1.5975 0.205284
## stimulus 34.3 12 15.356 2136 397.8339 < 2.2e-16 ***
## focus:stimulus 0.3 24 15.356 2136 1.4808 0.062289 .
## alcohol:stimulus 0.2 12 15.356 2136 1.8983 0.030276 *
## focus:alcohol:stimulus 0.1 24 15.356 2136 0.5317 0.969485
## stage 0.5 2 8.437 356 9.7672 7.421e-05 ***
## focus:stage 0.1 4 8.437 356 1.0657 0.373312
## alcohol:stage 0.0 2 8.437 356 0.2461 0.781995
## focus:alcohol:stage 0.0 4 8.437 356 0.0856 0.986846
## stimulus:stage 0.5 24 18.452 4272 4.8156 9.804e-14 ***
## focus:stimulus:stage 0.2 48 18.452 4272 0.8259 0.798308
## alcohol:stimulus:stage 0.2 24 18.452 4272 1.5525 0.041858 *
## focus:alcohol:stimulus:stage 0.2 48 18.452 4272 0.7592 0.887985
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## stimulus 0.00230 0.0000e+00
## focus:stimulus 0.00230 0.0000e+00
## alcohol:stimulus 0.00230 0.0000e+00
## focus:alcohol:stimulus 0.00230 0.0000e+00
## stage 0.58653 3.1174e-21
## focus:stage 0.58653 3.1174e-21
## alcohol:stage 0.58653 3.1174e-21
## focus:alcohol:stage 0.58653 3.1174e-21
## stimulus:stage 0.00023 0.0000e+00
## focus:stimulus:stage 0.00023 0.0000e+00
## alcohol:stimulus:stage 0.00023 0.0000e+00
## focus:alcohol:stimulus:stage 0.00023 0.0000e+00
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## stimulus 0.43240 < 2.2e-16 ***
## focus:stimulus 0.43240 0.1380588
## alcohol:stimulus 0.43240 0.0893189 .
## focus:alcohol:stimulus 0.43240 0.8740470
## stage 0.70748 0.0005176 ***
## focus:stage 0.70748 0.3619277
## alcohol:stage 0.70748 0.7031497
## focus:alcohol:stage 0.70748 0.9624416
## stimulus:stage 0.52543 3.705e-08 ***
## focus:stimulus:stage 0.52543 0.7126193
## alcohol:stimulus:stage 0.52543 0.0944288 .
## focus:alcohol:stimulus:stage 0.52543 0.7985395
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## stimulus 0.4469239 4.817097e-240
## focus:stimulus 0.4469239 1.351004e-01
## alcohol:stimulus 0.4469239 8.677605e-02
## focus:alcohol:stimulus 0.4469239 8.789647e-01
## stage 0.7114920 5.039157e-04
## focus:stage 0.7114920 3.621441e-01
## alcohol:stage 0.7114920 7.044778e-01
## focus:alcohol:stage 0.7114920 9.629908e-01
## stimulus:stage 0.5681699 1.156033e-08
## focus:stimulus:stage 0.5681699 7.226680e-01
## alcohol:stimulus:stage 0.5681699 8.748527e-02
## focus:alcohol:stimulus:stage 0.5681699 8.098047e-01# need to take effects apart and set contrasts
# the main effect of alcohol
alcohol_emm <- emmeans::emmeans(itt_afx, ~alcohol, model = "multivariate")
alcohol_emm # shows us the means## alcohol emmean SE df lower.CL upper.CL
## alcoholic 0.887 0.00777 178 0.872 0.903
## non-alcoholic 0.919 0.00783 178 0.903 0.934
##
## Results are averaged over the levels of: focus, stage, stimulus
## Confidence level used: 0.95# the main effect of stage
stage_emm <- emmeans::emmeans(itt_afx, ~stage, model = "multivariate")
stage_emm## stage emmean SE df lower.CL upper.CL
## stage.1 0.909 0.00632 178 0.896 0.921
## stage.2 0.908 0.00553 178 0.897 0.919
## stage.3 0.892 0.00637 178 0.879 0.904
##
## Results are averaged over the levels of: focus, alcohol, stimulus
## Confidence level used: 0.95emmeans::contrast(stage_emm, method = "trt.vs.ctrl", ref = 1, adjust = "holm")## contrast estimate SE df t.ratio p.value
## stage.2 - stage.1 -0.000635 0.00473 178 -0.134 0.8933
## stage.3 - stage.1 -0.017376 0.00543 178 -3.198 0.0033
##
## Results are averaged over the levels of: focus, alcohol, stimulus
## P value adjustment: holm method for 2 tests#main effect of stimulus
stimulus_emm <- emmeans::emmeans(itt_afx, ~stimulus, model = "multivariate")
stimulus_emm## stimulus emmean SE df lower.CL upper.CL
## stim1 0.701 0.00747 178 0.686 0.716
## stim2 0.846 0.00784 178 0.830 0.861
## stim3 0.857 0.00733 178 0.843 0.871
## stim4 0.862 0.00770 178 0.846 0.877
## stim5 0.921 0.00682 178 0.908 0.934
## stim6 0.924 0.00623 178 0.911 0.936
## stim7 0.934 0.00610 178 0.922 0.946
## stim8 0.944 0.00574 178 0.933 0.955
## stim9 0.945 0.00586 178 0.934 0.957
## stim10 0.948 0.00586 178 0.936 0.959
## stim11 0.947 0.00597 178 0.935 0.959
## stim12 0.954 0.00550 178 0.943 0.965
## stim13 0.956 0.00566 178 0.945 0.967
##
## Results are averaged over the levels of: focus, alcohol, stage
## Confidence level used: 0.95emmeans::contrast(stimulus_emm, method = "poly", adjust = "holm")## contrast estimate SE df t.ratio p.value
## linear 2.76 0.0932 178 29.566 <.0001
## quadratic -5.21 0.1856 178 -28.062 <.0001
## cubic 1.32 0.0952 178 13.873 <.0001
## quartic -7.12 0.9739 178 -7.307 <.0001
## degree 5 2.40 0.2484 178 9.660 <.0001
## degree 6 -3.96 0.3900 178 -10.149 <.0001
##
## Results are averaged over the levels of: focus, alcohol, stage
## P value adjustment: holm method for 6 tests# interaction stimulus:stage
inter_emm <- emmeans::emmeans(itt_afx, c("stage", "stimulus"), model = "multivariate")
inter_emm## stage stimulus emmean SE df lower.CL upper.CL
## stage.1 stim1 0.672 0.01041 178 0.652 0.693
## stage.2 stim1 0.716 0.01077 178 0.695 0.738
## stage.3 stim1 0.715 0.01037 178 0.694 0.735
## stage.1 stim2 0.863 0.00923 178 0.844 0.881
## stage.2 stim2 0.847 0.00921 178 0.829 0.865
## stage.3 stim2 0.828 0.00976 178 0.809 0.847
## stage.1 stim3 0.865 0.00876 178 0.848 0.883
## stage.2 stim3 0.870 0.00889 178 0.852 0.887
## stage.3 stim3 0.836 0.00896 178 0.818 0.854
## stage.1 stim4 0.881 0.00909 178 0.863 0.899
## stage.2 stim4 0.855 0.00907 178 0.837 0.873
## stage.3 stim4 0.849 0.00973 178 0.830 0.868
## stage.1 stim5 0.924 0.00835 178 0.907 0.940
## stage.2 stim5 0.928 0.00749 178 0.914 0.943
## stage.3 stim5 0.911 0.00841 178 0.894 0.927
## stage.1 stim6 0.936 0.00682 178 0.922 0.949
## stage.2 stim6 0.928 0.00727 178 0.914 0.943
## stage.3 stim6 0.907 0.00783 178 0.891 0.922
## stage.1 stim7 0.941 0.00780 178 0.925 0.956
## stage.2 stim7 0.937 0.00638 178 0.924 0.949
## stage.3 stim7 0.926 0.00749 178 0.911 0.940
## stage.1 stim8 0.951 0.00750 178 0.936 0.966
## stage.2 stim8 0.948 0.00575 178 0.937 0.959
## stage.3 stim8 0.932 0.00814 178 0.916 0.948
## stage.1 stim9 0.951 0.00758 178 0.936 0.966
## stage.2 stim9 0.953 0.00625 178 0.941 0.965
## stage.3 stim9 0.932 0.00689 178 0.919 0.946
## stage.1 stim10 0.956 0.00723 178 0.941 0.970
## stage.2 stim10 0.955 0.00653 178 0.942 0.968
## stage.3 stim10 0.932 0.00727 178 0.918 0.947
## stage.1 stim11 0.952 0.00783 178 0.937 0.968
## stage.2 stim11 0.953 0.00610 178 0.941 0.965
## stage.3 stim11 0.936 0.00765 178 0.921 0.952
## stage.1 stim12 0.963 0.00717 178 0.949 0.977
## stage.2 stim12 0.955 0.00566 178 0.944 0.966
## stage.3 stim12 0.943 0.00713 178 0.929 0.957
## stage.1 stim13 0.962 0.00665 178 0.949 0.975
## stage.2 stim13 0.962 0.00626 178 0.950 0.974
## stage.3 stim13 0.944 0.00715 178 0.930 0.958
##
## Results are averaged over the levels of: focus, alcohol
## Confidence level used: 0.95emmeans::contrast(
inter_emm,
c(stage = "trt.vs.ctrl", stimulus = "poly"),
rel = 1,
adjust = "holm"
)## contrast estimate SE df t.ratio p.value
## stage.2 stim1 - stage.1 stim1 0.0439 0.0135 178 3.257 0.0027
## stage.3 stim1 - stage.1 stim1 0.0421 0.0137 178 3.079 0.0027
## stage.1 stim2 - stage.1 stim1 0.1902 0.0106 178 18.003 <.0001
## stage.2 stim2 - stage.1 stim1 0.1742 0.0116 178 15.007 <.0001
## stage.3 stim2 - stage.1 stim1 0.1557 0.0117 178 13.329 <.0001
## stage.1 stim3 - stage.1 stim1 0.1930 0.0101 178 19.207 <.0001
## stage.2 stim3 - stage.1 stim1 0.1973 0.0104 178 18.880 <.0001
## stage.3 stim3 - stage.1 stim1 0.1634 0.0111 178 14.660 <.0001
## stage.1 stim4 - stage.1 stim1 0.2082 0.0110 178 18.840 <.0001
## stage.2 stim4 - stage.1 stim1 0.1829 0.0113 178 16.166 <.0001
## stage.3 stim4 - stage.1 stim1 0.1763 0.0119 178 14.801 <.0001
## stage.1 stim5 - stage.1 stim1 0.2515 0.0110 178 22.820 <.0001
## stage.2 stim5 - stage.1 stim1 0.2560 0.0112 178 22.949 <.0001
## stage.3 stim5 - stage.1 stim1 0.2382 0.0109 178 21.847 <.0001
## stage.1 stim6 - stage.1 stim1 0.2634 0.0106 178 24.888 <.0001
## stage.2 stim6 - stage.1 stim1 0.2560 0.0113 178 22.664 <.0001
## stage.3 stim6 - stage.1 stim1 0.2344 0.0114 178 20.518 <.0001
## stage.1 stim7 - stage.1 stim1 0.2682 0.0107 178 24.987 <.0001
## stage.2 stim7 - stage.1 stim1 0.2641 0.0105 178 25.037 <.0001
## stage.3 stim7 - stage.1 stim1 0.2531 0.0112 178 22.671 <.0001
## stage.1 stim8 - stage.1 stim1 0.2786 0.0108 178 25.780 <.0001
## stage.2 stim8 - stage.1 stim1 0.2756 0.0102 178 26.988 <.0001
## stage.3 stim8 - stage.1 stim1 0.2599 0.0116 178 22.321 <.0001
## stage.1 stim9 - stage.1 stim1 0.2784 0.0110 178 25.341 <.0001
## stage.2 stim9 - stage.1 stim1 0.2806 0.0111 178 25.329 <.0001
## stage.3 stim9 - stage.1 stim1 0.2598 0.0111 178 23.461 <.0001
## stage.1 stim10 - stage.1 stim1 0.2832 0.0110 178 25.756 <.0001
## stage.2 stim10 - stage.1 stim1 0.2827 0.0106 178 26.659 <.0001
## stage.3 stim10 - stage.1 stim1 0.2598 0.0117 178 22.153 <.0001
## stage.1 stim11 - stage.1 stim1 0.2798 0.0114 178 24.602 <.0001
## stage.2 stim11 - stage.1 stim1 0.2806 0.0108 178 25.919 <.0001
## stage.3 stim11 - stage.1 stim1 0.2640 0.0117 178 22.576 <.0001
## stage.1 stim12 - stage.1 stim1 0.2907 0.0109 178 26.560 <.0001
## stage.2 stim12 - stage.1 stim1 0.2828 0.0112 178 25.360 <.0001
## stage.3 stim12 - stage.1 stim1 0.2706 0.0117 178 23.037 <.0001
## stage.1 stim13 - stage.1 stim1 0.2892 0.0108 178 26.842 <.0001
## stage.2 stim13 - stage.1 stim1 0.2896 0.0109 178 26.655 <.0001
## stage.3 stim13 - stage.1 stim1 0.2714 0.0115 178 23.677 <.0001
##
## Results are averaged over the levels of: focus, alcohol
## P value adjustment: holm method for 38 testsitt_dat$stimulus <- as.factor(itt_dat$stimulus) %>%
fct_relevel(.,c("stim1", "stim2", "stim3", "stim4", "stim5", "stim6", "stim7", "stim8", "stim9", "stim10", "stim11", "stim12","stim13"))
itt_plot<-ggplot(data = itt_dat, aes(x = stimulus, y = correct, color = alcohol))+
facet_grid(rows = vars(focus),
cols = vars(stage))+
stat_summary(fun=mean, geom="point") +
scale_y_continuous(name="% correct", limits=c(0.5, 1.0), breaks = seq(0.5,1,0.15)) +
scale_x_discrete(name = "stimulus duration\n(ms)", labels =c("19", "25", "31", "37", "44", "50", "62", "75", "87", "100", "125", "150", "200") )+
scale_color_manual(values = gender_pal) +
theme(axis.text.x = element_text( color="black",
size=2, angle=45))+
theme_minimal()
ggplotly(itt_plot, height = 500, width=1200)Participants’ cognitive performance
mood_data<- mood_dat %>%
dplyr::select(subject, stage, calm, content, alert)%>%
pivot_longer(cols = c(calm, content, alert),
names_to = "factor",
values_to = "rating")%>%
mutate(condition = subject%/% 1000)%>%
mutate(alcohol = if_else((condition%%2)==1, "alcoholic", "non-alcoholic"))%>%
mutate(focus = case_when(condition == 1 ~ 'control',
condition == 2 ~ 'control',
condition == 3 ~ 'hedonic',
condition == 4 ~ 'hedonic',
condition == 5 ~ 'utilitarian',
condition == 6 ~ 'utilitarian'))
ggplot(data = mood_data, aes(x = stage, y = rating, color = alcohol))+
facet_grid(rows = vars(focus),
cols = vars(factor))+
stat_summary(fun=mean, geom="point") +
scale_y_continuous(name="rating", limits=c(35, 48), breaks = seq(0,50,5)) +
scale_x_discrete(name = "stage", labels = c("1", "2", "3") )+
scale_color_manual(values = gender_pal) +
theme(axis.text.x = element_text( color="black",
size=2, angle=45))+
theme_light()
Participants’mood accross the stages and conditions
Assessing the effects of alcohol content and focus on mood, we observed main effect of mood type, main effect of stage. We have also observed interactions between focus and alcohol, mood and stage and a three-way intearction between mood, alcohol and stage. We will briefly report all main effects and interactions but will focus on the interpretation of the three-way interaction and the two-way intercation between focus and alcohol content.
There was a significant main effect of factor (from the Bond-Lader questionnaire: calm, content, alert) on ratings. Participants’ ratings of feeling calm were significantly lower than feelings related to content and alertness. Participants also rated feelings of alertness significantly higher than feelings of calm or content. Feeling content was rated significantly higher than feeling calm, but significantly lower than feeling alert. careful–> interactions!!
There was a significant effect of stage on mood ratings. Mood (all three factors) was rated significantly higher in stage 3 than in stage 1, there was no difference between stage 1 and stage 2 or stage 2 and stage 3. careful–> interaction!!
The ratings were affected by focus and alcohol content of the beer. In the hedonic condition, compared to control condition, the ratings were significantly higher for participants who drank alcoholic beer. This difference, however was not present in the utilitarian condition.
Interestingly, participants self-reported feeling more alert at the end (stage 3) of the experiment than at the beginning (stage 1). Participants also rated feeling less calm during the experiment (stages 2 and 3) than at the beginning. This was regardless of the conditions they were in. careful interaction
Participants rated feeling more alert as the experiment progressed, more so in the alcoholic beer condition.